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Abstract 


A new algorithm is presented to integrate component balances along polymer electrolyte membrane fuel cell (PEMFC) channels to obtain 
three-dimensional results from a detailed two-dimensional finite element model. The analysis studies the cell performance at various hydrogen 
flow rates, air flow rates and humidification levels. This analysis shows that hydrogen and air flow rates and their relative humidity are critical 
to current density, membrane dry-out, and electrode flooding. Uniform current densities along the channels are known to be critical for thermal 
management and fuel cell life. This approach, of integrating a detailed two-dimensional across-the-channel model, is a promising method for fuel 
cell design due to its low computational cost compared to three-dimensional computational fluid dynamics models, its applicability to a wide range 
of fuel cell designs, and its ease of extending to fuel cell stack models. 


© 2006 Elsevier B.V. All rights reserved. 
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1. Introduction 


Fuel cell design and technology are critical research areas 
necessary for the transition to a hydrogen economy. There 
are several different fuel cell technologies, polymer electrolyte 
membrane fuel cell (PEMFC), alkaline fuel cell (AFC), phos- 
phoric acid fuel cell (PAFC), molten carbonate fuel cell (MCFC) 
and solid oxide fuel cell (SOFC) each with different designs, 
chemistries and operating conditions [1]. PEMFC’s stand out 
from other fuel cell technologies due to their higher current 
densities, use of a solid polymer membrane for the cell elec- 
trolyte, and low temperature operation (60-90°C). By using 
a solid polymer electrolyte the problems associated with seal- 
ing, handling, and assembly of corrosive electrolytes used in 
the other fuel cell technologies are avoided. Furthermore the 
low temperature operation of PEMFC makes them suitable for 
use in mobile, transportation and backup power systems due to 
their fast start-up times, easier thermal management and overall 
system weight [1,2]. The major challenges facing commercial- 
ization of the PEMFC are the cost, reliability and life [1,3]. 
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Mathematical models and simulations are used extensively 
by researchers to overcome the engineering challenges of com- 
mercializing PEMFC’s. A significant amount of research has 
been aimed at increasing the understanding of PEMFC physics 
through modeling. It is beyond the scope of this paper to review 
all of the models; a detailed review of the available fuel cell mod- 
els and their comparisons is available elsewhere (Haraldsson and 
Wipke [4]). 

The models developed by researchers have evolved from one- 
dimensional “along-the-channel” forms [5-7] to more realistic 
two-dimensional models of the flow channels and a fundamen- 
tal representation of the membrane electrode assembly (MEA) 
(Gurau et al. [8]). With the aid of computational fluid dynamics 
(CFD) software packages in late 1990s, the fuel cell modelers 
began developing three-dimensional CFD models [9-13]. The 
earlier two-dimensional models were mostly along-the-channel 
models and the effects of bipolar plate shoulders and chan- 
nel sizes were not studied. The three-dimensional CFD models 
solved this limitation by adding the across the channel third 
dimension. 

Commercial fuel cell anode and cathode catalyst layer thick- 
nesses are on the order of 10-20 wm and the height and width 
of the single cell in a stack are 10-20 cm. These dimension 
mismatches require at least a million elements to model even 
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Nomenclature 

Cw mass concentration of water in the membrane 
(kg m7) 

Dy water diffusivity (m? s7!) 

Fi flow rate 

F Faraday’s constant, 96487 C (mol~!) 

na drag coefficient 

Nw water flux (mol m~? s7!) 

p pressure (Pa) 

t thickness (m) 

Veell cell operating potential (V) 

W width (m) 

x mole fraction 


Greek letters 


o conductivity (S m7!) 

(0 potential (V) 

E porosity 

Subscripts 

a anode 

c cathode 

GDE gas distribution electrode 

i components, H2, and H20 for the anode, O2, H20, 


and N3 for the cathode 
liquid water 


— 


m membrane 

w water in the membrane 
Superscripts 

0 boundary condition 


a small 7cm x 1 cm section of a fuel cell, and can take more 
than 1 h with a 10-50 node parallel computer [13]. Although the 
recent modeling trend is three-dimensional CFD PEMFC mod- 
els, in our previous work [14], we developed a two-dimensional 
across-the-channel CFD PEM fuel cell model where the effects 
of channel and bipolar plate dimensions as well as the properties 
of the fuel cell components were studied for various given con- 
ditions. By using a two-dimensional across-the-channel model 
we were able to study the effects of cell design and operating 
conditions rapidly and accurately. 

In later work [15], we used the two-dimensional finite element 
model and a design-of-experiments approach to study the main 
and interaction effects of five key design factors. These factors 
included: channel size, shoulder size, gas distribution electrode 
(GDE) thickness, GDE porosity, and GDE conductivity. These 
were all studied at moderate and high current density operations. 
An interaction is defined as the failure of a factor to produce the 
same response at different levels of another factor. Strong inter- 
actions between design factors that effect mass transport were 
identified and quantified which will aid researchers in the design 
of PEMFCs. Because the approach used was an across-the- 
channel model, it did not study the concentration change effects 


along the channel. To observe those effects the analysis was run 
for the inlet and exit fuel cell gas concentrations. The results of 
that work showed that the interactions between design factors 
become more significant at the exit conditions of the fuel cell. 

For a given design, fuel cell performance depends strongly on 
operating conditions. The effects of operating temperature and 
pressure have been studied by various researchers with mathe- 
matical models [9,16]. The range of these operating parameters 
is limited by the physical and chemical limitations of the PEMFC 
materials. 

To achieve the maximum power thus reducing the cost per 
kW, a major challenge is water management. The fuel cell must 
be operated such that the membrane is adequately hydrated to 
reduce the ohmic losses. Also the water generated at the cathode 
must be removed quickly too avoid electrode flooding. Recent 
work by others has focused on humidification levels and flow 
rates of hydrogen and air. Berning and Djilali [17] developed a 
two-phase (liquid—gas) model, of the porous GDE and the gas 
channel of a PEMFC but excluded the MEA. Water transport 
through the MEA was therefore ignored. Recently, Pasaogullari 
and Wang [18] conducted a comprehensive study on the stoi- 
chiometry and humidification requirements of PEMFC opera- 
tion using a two-phase model. Their work showed that flow rates 
and humidification are key parameters for highly efficient fuel 
cells. The ratio of reactant available divided by the amount of 
reactant required for the desired current density (stoichiometric 
ratio), is also a critical operating condition as it determines the 
efficiency of the fuel cell. Feed flow rates of hydrogen and water 
are also critical parameters to manage water in the fuel cell itself 
and fuel cell system [19]. Berning and Djilali [17] also studied 
the effects of stoichiometric ratio on performance but because 
their model assumes constant membrane water content, it could 
not answer the water management questions. 

Even though our two-dimensional across-the-channel finite 
element model [14] is fast and detailed, it lacks the ability to 
study the effects of concentration changes along the channel. 
This is a significant limitation when studying the flow rate and 
humidification requirements. In this study, a simple Euler inte- 
gration method has been added to the two-dimensional model to 
solve the problem of flow of reactants along the channel. With 
the algorithm developed in this study the performance of the 
fuel cell can be evaluated in a fraction of the time of a compa- 
rable three-dimensional model. This increased speed allows the 
model to be used as a design tool, not just an analysis tool. 


2. Two-dimensional model 


The across-the-channel model used in this study was pub- 
lished earlier [14]. It assumes: 


Steady-state operation, 

isothermal operation, 

ideal gas mixtures, 

single phase model, 

isotropic and homogeneous electrodes and membrane, 

the membrane is considered impermeable for the gas phase, 
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channel exit, and channel-GDE boundary. Control volumes for 
both anode and cathodes are chosen that include the channels 
and GDE’s at each side of the channel. 

The algorithm for solving the problem is shown in Fig. 2. To 
simulate a three-dimensional model first the inlet flow rate of 
components are calculated. Because the velocity profiles in the 
channels are not modeled, the pressure drop along-the-channel is 
neglected. Due to the high flow rates and small channels widths 
and heights, concentration gradients across the gas channel (x- 
and y-directions) are neglected. The transport of the components 
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0 0 
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Fig. 1. The dependent variables, boundary conditions, and computational 
domain. 


e negligible contact resistance, 
e minimal membrane swelling. 


The anode and cathode catalyst layers are modeled as reactive 
boundaries because their thickness of 10-25 um is significantly 
smaller than all other component thicknesses. The remainder of 
the fuel cell model is a comprehensive two-dimensional, isother- 
mal, steady-state model providing a detailed description of the 
following transport phenomena: 


Multi-component flow, 

diffusion of reactants through the porous electrodes, 
electrochemical reactions, 

transport of electrons through the electrodes, 

water balance in the membrane. 


The equations governing these processes include: 


e Ionic balance in electrodes and membrane, Gpe and ¢m, 

e the Maxwell—Stefan equations for multi-component diffusion 
and convection in gas distribution channels and gas distribu- 
tion layers, xH,, xo,, and xH,0,c, 

e Darcy’s law for the flow of species in porous electrodes, p, 

e water balance and water flux in the membrane governed by 
diffusion, convection and electro-osmotic drag, pı and cy. 


The domain, dependent variables and boundary conditions 
are outlined in Fig. 1, details of the model can be found in Guve- 
lioglu and Stenger [14]. 

FEMLAB® a finite element computational fluid dynamics 
package, was used to solve the nonlinear system of equations. 
The model was first created and tested in FEMLAB®’s graph- 
ical user interface and then saved as a MATLAB® m-file [20]. 
Details of the model, geometry meshing and solver settings were 
presented previously [14]. 


3. Two-dimension to three-dimension integration 


The anode and cathode channels have three face boundaries 
where components can enter or leave, these are channel entrance, 


Input 
Fuel and Oxidant 
e Composition 
e Relative Humidity 
e = Flow Rates 


Calculate Inlet Molar Flow Rates 
Anode Cathode 


Ny, (0),Ni,0(9) No, (0), Nu,o(0), Ny, (0) 


dN y, (K) dNy,o(k) 


Calculate Mole Fractions 
Anode Cathode 


XH, >XH 0o Yo, »¥m,o> N, 


Call FEMLAB® 2D Model in 
MATLAB® 


Calculate 7”, on s 


Anode Cathode 
ANo, (k) dN yolk) 


dz dz dz dz 


Update Molar Flow Rates by Simple 
Euler Integration 


Nerd- e 
dz 


z(k+1)=z(k)+h 


Fig. 2. Two-dimension to three-dimension integration algorithm. 


G.H. Guvelioglu, H.G. Stenger / Journal of Power Sources 163 (2007) 882-891 885 


in the GDE’s as well as water in the membrane along the channel 
(z-direction) is neglected due to the order of magnitude differ- 
ence in dimensions of the channel width and height versus the 
length. 

The inlet mass flow rates are calculated based on the specified 
flow rate, where F, is the anode flow rate and is equal to 1.0 when 
the hydrogen flow rate is sufficient to produce 1 Acm7?. 

For the anode, the inlet flow rate of hydrogen is calculated 
with the following equation: 


1 x 104 x (Wen + Wen) x L 


Nu, (0) = Fa x IST a) 


where W,, and Wen are the shoulder and channel widths, L the 
length of the channel and F is Faraday’s constant. 

For inlet compositions, the water mole fraction is calculated 
from the relative humidity (RH): 

sat X RH, 
tmoa = AA (2) 
Pa 

where RH, is the relative humidity of the anode stream, pa the 
anode operating pressure, Psat is the saturation pressure of water 
calculated with [21]: 


7258.2 
Psat = exp (73.688 a 7 7.3037 logT 


+ 4.1653 x 10-°T? (3) 


After the water mole fraction at the anode is calculated the molar 
flow rate of water is calculated by Eq. (4): 


xmo X Np, 


N (0) = 
H50,a(0) ae 


(4) 


The cathode flow rates are calculated similar to the anode. 


The flow rate, Fo is used and is equal to 1.0 when the oxygen 


flow rate is sufficient to generate 1 Acm™~°. 


1 x 10* x (Wa + W, L 
No,(0) = F; x x x ` sh + Wen) X 5) 
xF 


The inlet nitrogen molar flow rate is calculated using oxygen 
molar flow rate and the composition of air before humidification. 


No, x 0.79 
Nn, (0) = ————— 6 
No ( ) 0.21 ( ) 
The water mole fraction is calculated by: 
sat X RH 
amo = A= (7) 
Pc 


And the inlet water molar flow rate for the cathode side is 
calculated by: 


XH20,c X (No, + Nn») 
(1 — xy,0) 


Nn,o,c(0) = (8) 

Once the flow rates of the components are calculated the mole 
fractions of components are used as boundary conditions for the 
GDE-gas channel interface for the two-dimensional model as 


shown in Fig. 1. The two-dimensional model is then used to 
solve mass, momentum and charge balances. The eight depen- 
dent variables, p, pi, Cw, XH), XO, XH20,c> PGDE and dy are then 
obtained in the computational domain. 

A control volume is selected that includes both the gas chan- 
nel and the GDE. The total flux of hydrogen and water at the 
anode catalyst-GDE boundary and the total flux of oxygen and 
water at the cathode catalyst-GDE boundary (Ws) + Wen) as 
shown in Fig. 1, is calculated by integrating the flux values along 
the boundary. 


dN; k Wont Wen 
“A =| Ji(y) dy (9) 


where j; is the molar flux of component i, and i is equal to H and 
H20 for the anode and O2 and H20 for the cathode. The gradient 
of the flux of N2 is zero because the membrane is impermeable 
to all components except water. The average current density 
is calculated by integrating the local current density along the 
catalyst layer as shown in Eq. (10). 


1 Weh+Weh 
pret | 1()) dy (10) 
Wen + Wen Jo 


For further insight, the average conductivity of the mem- 
brane is calculated by integrating the membrane conductivity 
and dividing it by the membrane area. 


1 Wen+Wen tm 
ont = a f Om(x, y) - dx - dy 
i (Wh T Wen) ‘tm Jo 0 E 


(11) 
Table 1 
Base case geometric parameters 
Parameter Value 
Wen, the gas channel width (m) 7.62 x 1074 
Wsn, the bipolar shoulder width (m) 7.62 x 1074 
L, channel length (m) 0.1 
tape, anode and cathode GDE thickness (m) 3 x 1074 
tm, membrane thickness (m) 1.78 x 1074 
€GpE, GDE porosity of the anode and cathode 0.6 
Table 2 
Operating condition parameters 
Symbol Value 
Veen, cell potential 0.6 V 
T, temperature 343.15K 
pe, anode side pressure 202650 Pa 
pe, cathode side pressure 202650 Pa 

Relative humidity 
100% 75% 50% 

xe, , anode feed hydrogen mole fraction 0.846 0.885 0.923 
Seas anode feed water mole fraction 0.154 0.115 0.077 
PA , cathode feed oxygen mole fraction 0.178 0.186 0.194 
Woa cathode feed water mole fraction 0.154 0.115 0.077 
me , cathode feed nitrogen mole fraction 0.668 0.699 0.729 


886 G.H. Guvelioglu, H.G. Stenger / Journal of Power Sources 163 (2007) 882-891 


The fluxes of the components in and out of the GDE-channel 
control volume and the inlet mass flow rates create a boundary 
value problem in the third dimension (z). A simple variable step 
size Euler integration algorithm is used to solve this problem by 
integrating along the channel in z-direction. 


dNj(k) 


Nik + 1) = Ni(k) +h x (12) 
where kis the grid point and i = H2 and H20 for the anode and O2 
and H20 for the cathode and A is the step size in the z-direction. 

The step size along the channel, h is determined such that 
the relative change of the molar flow rates of components are 
less than 2%. After step size refinements and accuracy checks, a 
2% relative change criteria gave a reasonable balance between 
accuracy and computing time. 

In addition to concentration and current profiles along the 
channel, overall power output is calculated by integrating the 
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~ 


current density along the channel using the trapezoid rule and 
multiplying by the cell potential. 


k-1 


h x (IYS(k) + IYS(k + 1) 


5 (13) 


Power = Veell X 


k=0 
Overall power output is then used to calculate power density 
by dividing by the effective MEA area: 


power 


(Wen + Wen) x L oe 


power density = 

This integration algorithm is implemented in MATLAB® 
as an M-File which calls the FEMLAB® two-dimensional 
model. The integration along the channel required 30-70 two- 
dimensional model calls. The typical run times for each 2D 
simulation are 30-100 s, with the longer times for high current 
density operation. The total simulation times range between 30 
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Fig. 3. Profiles along the channel for 100% relative humidity fuel and oxidant: (a) average current density, (b) average anode hydrogen mole fraction, (c) average 
cathode oxygen mole fraction, (d) average cathode water mole fraction, (e) average membrane conductivity, (f) average net water flux density. 
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and 70min. The majority of the simulations for this problem 
were completed under 45 min on an Intel Pentium® 4 3.2 Ghz 
CPU with 1 GB of DDRam. 


4. Results and discussions 


The model was run for flow rates of 1 and 3 based at 1 Acm7?, 
for both the anode (F,) and the cathode sides (Fe), at three differ- 
ent humidification levels (100, 75, and 50% relative humidity) 
and at a fuel cell potential of 0.6 V. The geometric parameters of 
the fuel cell used are listed in Table 1. The fuel and oxidant flows 
are co-current in the flow channels. The operating conditions of 
the cell and the inlet composition of the fuel and oxidant are 
shown in Table 2. Fig. 3 shows the results for 100% RH hydro- 
gen and air, Fig. 4 shows the results for 75% RH hydrogen and 
air, and Fig. 5 has the results for 50% RH case. 
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The average current density profile, calculated by Eq. (10), 
along a 10cm long straight channel is shown in Fig. 3(a) for 
four different stoichiometric flow rates. This plot shows that 
increasing the hydrogen flow rate (F,) from 1 to 3 increases 
the current density significantly, mostly because increasing the 
hydrogen flow rate also increases the amount of water supplied. 
More water in the anode feed hydrates the membrane causing 
higher membrane conductivity; Fig. 3(e) shows the increased 
membrane conductivity. 

Fig. 3(b) is a plot of hydrogen mole fraction versus length. 
This confirms that the water concentration does increase with 
hydrogen flow (Fa =3). 

Increasing the air flow rate (Fs = 3) does not improve the cur- 
rent density but can actually decrease the performance slightly 
as shown in Fig. 3(a). Even though higher air flow maintains 
a higher oxygen concentration, the higher air flow prevents the 
concentration of water to increase on the cathode side. This 
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Fig. 4. Profiles along the channel for 75% relative humidity fuel and oxidant: (a) average current density, (b) average anode hydrogen mole fraction, (c) average 
cathode oxygen mole fraction, (d) average cathode water mole fraction, (e) average membrane conductivity, (f) average net water flux density. 
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Fig. 5. Profiles along the channel for 50% relative humidity fuel and oxidant: (a) average current density, (b) average anode hydrogen mole fraction, (c) average 
cathode oxygen mole fraction, (d) average cathode water mole fraction, (e) average membrane conductivity, (f) average net water flux density. 


is caused by the air removing water from the cathode chan- 
nel (Fig. 3(d)). Lower water concentrations at the cathode side, 
causes less back diffusion of water from cathode to anode and 
thus lowers the membrane conductivity. 

When the PEM fuel cell is operated at equal pressures at 
anode and cathode sides, water transport in the membrane is 
governed by only two processes: flow of water from anode to 
cathode due to electro-osmotic drag and the diffusion of water 
from higher to lower concentrations. The convective flux of 
water is zero inside the membrane due to the zero pressure gra- 
dient. Thus the water flux can be written as: 


na- I 


Nw = — DwVcw (15) 


where cw is the concentration of water, Dw the diffusion coef- 
ficient of water in the membrane, ng the electro-osmotic drag 


coefficient, I the local current density vector, and F is the Fara- 
day constant. 

Fig. 3(b) shows that the anode water mole fractions are 
between 0.09 and 0.15 (1 — anode hydrogen mole fraction), 
however on the cathode side the mole fraction of water ranges 
between 0.154 and 0.375 (Fig. 3 (d)). The net water flux in 
Fig. 3(f) is always positive meaning flow of water is from the 
anode to the cathode. 

When the current density profiles for all four cases in Fig. 3(a) 
are compared it is interesting that the current densities converge 
to almost the same value at the exit. However hydrogen mole 
fraction values converge to a single value at the exit only for a 
given air flow rate (Fig. 3(b)). For example at Fo = 1 the hydro- 
gen mole fractions converge to 0.893, and at Fo =3 the hydrogen 
mole fractions converge to 0.906. This shows that water concen- 
trations at the anode and the cathode sides reach equilibrium and 
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the equilibrium value appears to be a slight function of the air 
flow rates. 

A lower equilibrium water concentration at the anode for the 
higher air flow rates S, = 3 also confirms the water removal effect 
of air. Fig. 3(c) shows the mole fraction profile of oxygen along 
the cathode channel. The mole fraction of oxygen is almost linear 
along the channel. As expected, higher air flow rates result in 
a slower decrease of oxygen and higher exiting mole fractions. 
The oxygen mole fraction profiles show slight differences for 
the two hydrogen flow rates. The higher hydrogen flow rates 
(Fa =3) causes more water to be transported to the cathode side 
thus lowering the mole fraction of oxygen. 

The 75% RH hydrogen and air results are shown in Fig. 4. 
Similar to the 100% RH case, the cell performance can be 
improved by increasing the hydrogen flow rate as shown in 
Fig. 4(a). Also the water removal effect of the higher air flow 
rates on the cell performance is more noticeable for 75% RH 
(Fig. 4(d)) than for the 100% RH case due to lower water supply. 

Fig. 4(a) shows that the current density at the entrance for the 
75% RH cases is significantly lower than the 100% RH cases 
(Fig. 3(a)). This is caused by the higher membrane conductivity 
at the entrance 10.3S m7! for 100% RH versus 8 S m7! for the 
75% relative humidity (Figs. 3(e) and 4(e)). 

Similar to the 100% RH case, the 75% RH case reaches an 
equilibrium water flux at the anode exit. As expected for all flow 
rates the conductivities are always lower for all lower RH cases. 
For Fe =3 the back diffusion of water from cathode to anode is 
not high enough to keep the membrane hydrated at the entrance 
levels, due to the water removal effect of the higher air flow 
rate. For the F, = 1 cases, membrane conductivities are higher 
at the exit than the entrance due to higher water concentrations 
obtained at the exit. This increase causes an increase in current 
density at approximately 2cm from the entrance (Fig. 4(a)). 
Although the membrane conductivity starts increasing after the 
first 2 cm, it does not improve the current density significantly, 
because the concentration of oxygen also decreases rapidly at 
the exit. 

Fig. 4(f) shows the net water flux through the membrane for 
75% RH case. The net water flux at the entrance is less than half 
of the 100% RH case mainly because higher current densities 
at 100% RH, moves more water from anode to cathode due to 
electro-osmotic drag. 

Although the 100 and 75% RH case performances show simi- 
larities, the results for 50% RH are significantly different (Fig. 5). 
The entrance current density is 0.3 A cm~? for the 50% RH case 
which is significantly lower than the 100 and 75% RH cases 
of 0.70 and 0.58 Acm~”, respectively. This is due to the very 
low conductivity of the membrane at the entrance (2.9 S m-!). 
As the water produced at the cathode side increases, the mem- 
brane becomes more hydrated and its conductivity increases. 
This increases the current density along the channel as shown 
in Fig. 5(a). Similar to the 100 and 75% RH cases, increasing 
hydrogen flow rates improved the performance for the 50% RH 
cases for the first half of the channel; however after 5cm, the 
performance of Fa = 1 (Fe = 1) is better than Fa =3 (Fe = 1). For 
the 50% RH case, high air flow rates significantly decrease the 
current density and membrane conductivity due to the water 


removal effect of the high flow rates. In general the low over- 
all performance of 50% RH case is caused by the membrane 
not being properly hydrated, making it an undesired operating 
condition. 

For both the 100 and 75% RH cases the net water fluxes are 
from anode to cathode; however Fig. 5(f) shows the water flux 
changes from anode to cathode (positive) to cathode to anode 
(negative) after the first 3 cm of the channel. This is caused by the 
concentration of water increasing at the anode due to increased 
consumption of hydrogen causing back diffusion of water from 
the cathode. 


4.1. Flooding concerns 


At 343.15 K and 2 atm, water will condense when the mole 
fraction of water exceeds 0.154. Fig. 3(d) shows that water could 
condense throughout the cathode for 100% RH and it will start 
condensing after the first 1.3 cm of the fuel cell cathode channel 
for the 75% RH case. Pasaogullari and Wang [18] showed that 
electrode flooding will occur when current densities are higher 
than 1.4 A cm7? at 100% RH hydrogen and air. This is signif- 
icantly higher than the current densities studied in this work. 
Also in the Pasaogullari and Wang [18] work water activity lev- 
els (pw/p%') exceeded 3 in severe flooding cases, whereas the 
highest water activity level in this work is 2.4. It is uncertain 
that whether electrode flooding will occur in the cases shown in 
Figs. 3 and 4, however our higher GDE porosity of 0.6 versus 
0.5 of Pasaogullari and Wang [18], as well as our lower current 
densities, suggest that flooding may not occur. 

Figs. 6 and 7 show the power density calculated by Eq. (14), 
versus hydrogen utilization (varying F4) for 100 and 75% RH. 
For both 100 and 75% RH cases the power density decreases 
with increasing hydrogen utilization (decreasing Fa). The rel- 
ative magnitude of this decrease is higher for 100% RH than 
for 75% RH. Higher hydrogen flow rates (increasing Fa) show 
improved performance by providing more water and maintaining 
higher membrane hydration levels along the channel. However 
higher Fa’s decrease hydrogen utilization, requiring hydrogen 


3.45 I I 


3.4 — 


a 

oo 

a 
I 


iad 
w 
I 


o 

N 

a 
I 


3.2 


Power Density [kW m°] 


10 30 50 70 90 
Hydrogen Utilization [%] 


Fig. 6. Power density and hydrogen utilization for 100% RH fuel and oxidant. 
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Fig. 7. Power density and hydrogen utilization for 75% RH fuel and oxidant. 


recycling to maintain the fuel cell system efficiency. For both 100 
and 75% RH, increasing air flow rate decreases the power density 
more for high hydrogen utilizations. The effect of increasing air 
flow rates is less at low hydrogen utilizations because the power 
density is higher thus more water is generated at the cathode 
compared to high hydrogen utilization. With more water gen- 
erated at low hydrogen utilization the water removal effect of 
increasing air flow rates is smaller. 


5. Conclusions 


A new algorithm is presented to obtain three-dimensional 
results from a two-dimensional PEMFC model. A simple Euler 
integration method is used to solve the component balances 
along the fuel cell channels that uses a two-dimensional, across- 
the-channel, finite element model for local current density and 
water flux values. 

It was found that high hydrogen flow rates and low air flow 
rates performed better for both 100 and 75% RH cases. This is 
due to the water concentration in the membrane and the result- 
ing high membrane conductivity. For these RH cases the net 
water flux was from anode to cathode as the electro-osmotic 
drag exceeded the diffusive flux. Increasing the flow of hydro- 
gen improved the performance due to more water supplied to 
the anode side where it is needed most for humidification of the 
membrane. Higher hydrogen flow rates also resulted in lower uti- 
lization of hydrogen thus requiring hydrogen recycle for overall 
system efficiency. Increasing air flow rates lowered the perfor- 
mance for the studied conditions, as high air flow rates remove 
the generated water quickly and does not allow the concentration 
of water to build up at the cathode to hydrate the membrane. 

The low performance for hydrogen and air at 50% RH showed 
that the membrane is not properly hydrated. Even though the 
water concentration near the exit and the membrane conductivity 
increases, the performance of hydrogen and air with 50% RH is 
poor compared to 100 and 75% RH. 

The 100% RH case had a power density of 3.23 kW m~? at 
F,=1 with 50% hydrogen utilization; however, the maximum 


power density of 75% RH operation at Fs = 1 with 50% hydrogen 
utilization is 3.05 kW m7”. This is a 5.6% decrease in the power 
density due to lower humidification. As the hydrogen utilization 
increases, the difference between 100 and 75% RH operation 
power densities decreases. Operation at 75% RH gave a more 
uniform current density profile compared to 100% RH. Uniform 
current density in the PEMFC’s is desired to enhance thermal 
management, reliability and fuel cell life. Non-uniformities in 
current density can create hot spots due to local heat generation 
which can damage the membrane permanently. The 75% RH 
results show the trade off between performance and reliability 
which must be optimized to meet application requirements. 

The model used in these calculations assumes a single vapor 
phase thus liquid water formation and electrode flooding cannot 
be addressed. However when compared to the two-phase models 
reported in the literature that show electrode flooding, the water 
activity levels reached in this study are below reported flooding 
values. 

This study shows that the relative humidity of the hydrogen 
and the air, as well as their flow rates, are critical for the relia- 
bility, i.e. uniform current density, and power output. Air flow 
rate and relative humidity of air must be controlled to prevent 
membrane dry-out and electrode flooding. To study the extreme 
cases of electrode flooding a two-phase model should be used 
in the electrodes. 

The integration algorithm used in this work is the first attempt 
to combine a state-of-the-art finite element model with Euler 
integration to simulate a three-dimensional fuel cell. This new 
approach when compared to other three-dimensional CFD mod- 
els is much faster. With a detailed two-dimensional across-the- 
channel model, as used in this study, the effects of channel and 
shoulder dimensions can be studied. The order of magnitude 
differences in the characteristic dimensions of a PEMFC makes 
this approach more suitable. It has been shown in our previous 
work [14] that the concentration of reactants can change more 
than 90% under the bipolar plate shoulders for 1 mm channel 
and shoulder sizes. However due to the high flow rates in chan- 
nels only a 2% concentration change can be seen in 0.2-3 mm of 
channel length depending on operating conditions. Thus the rate 
of change of composition along the channel is small compared 
to that observed across the channel, validating our approach. 

This approach can be used for estimating overall PEMFC 
performance for co-current flow and parallel channel design. 
Also with an iterative procedure it could be used for counter- 
current flow as well as different channel layouts. 
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